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We report first-principles calculations for a ferroelectric organic crystal of phenazine and chlo- 
ranilic acid molecules. Weak intermolecular interactions are properly treated by using a second 
version of van der Waals density functional known as vdW-DF2 [K. Lee et a/., Phys. Rev. B 82, 
081101 (2010)]. Lattice constants, total energies, spontaneous electric polarizations, phonon modes 
and frequencies, and the energy barrier of proton transfer are calculated and compared with PBE 
and experiments whenever possible. We show that the donation of one proton from a chloranilic 
acid molecule to a neighboring phenazine molecule is energetically favorable. This proton transfer is 
the key structural change that breaks the centrosymmetry and leads to the ferroelectric structure. 
However, there is no unstable phonon associated with the proton transfer, and an energy barrier of 
8 meV is found between the paraelectric and ferroelectric states. 

PACS numbers: 61.66.Hq, 71.15.Mb, 77.80.-e, 78.55.Kz 



I. INTRODUCTION 



ization, estimated experimentally at 1-2/iC/cm 2 , is par- 



Ferroelectric materials are an important class of mate- 
rials. Their responses to various external stimuli can be 
used for many applications such as memory devices, elec- 
tromechanical actuators, ultrasonic sensors, electro-optic 
devices, and infrared thermal image sensors. Although 
transition-metal oxides are most widely used for applica- 
tions, organic ferroelectrics could be attractive alterna- 
tives, because they could be non- toxic, flexible, and easy 
to process i However, ferroelectric organic materials are 
rare, and substantial efforts are being made to find such 
materials that could be of practical use^ The cocrystal 
of phenazine (Phz) and chloranilic acid (H^ca) is one of 
several recently discovered hydrogen-bonded organic fer- 
roelectrics that have superior crystallinity and properties 
compared to conventional ferroelectric polymers^ - — 

The crystal structure of Phz-E^ca has been deter- 
mined by X-ray^i and neutron^ diffraction experiments. 
The centrosymmetric paraelectric structure (monoclinic 
P2i/ra, T > T c = 253 K) is shown in Fig. [TJ There are 
two molecules of each type (72 atoms in total) per unit 
cell. The two constituent molecules, Phz and Hbca, can 
make hydrogen bonds (H-bonds) with each other, form- 
ing linear chains in the crystal. One such chain that 
runs along the [110] direction in an ab plane is shown in 
Fig. DJb) as viewed along the b and c axes (top and bot- 
tom subpanels respectively). These chains stack along 
the b direction and fill an ab plane, as shown in Fig.[TJc). 
The next plane above or below [Fig. DJd)] is, however, 
filled with chains that run along [110]. 

Below 253 K, Phz-H^ca becomes ferroelectric^ In this 
polar phase (monoclinic P2i), one of the O-H bonds 
in H2ca stretches (by about 0.37 A) toward the N atom 
in the H-bonded Phz neighbor, adopting the structure 
shown in Fig.[2](b). 8 (The large arrows in this figure rep- 
resent the directions of the proton displacements. The 
other panels will be discussed later in Sec.lIIIi) The polar- 
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FIG. 1. (Color online) (a) Monoclinic unit cell of the paraelec- 
tric structure containing two H2ca molecules (at body-center 
and corner sites) and two Phz molecules (at a6-face-center 
and c-axis edge-center sites), (b) H-bonded chain of molecules 
running along [110]. (c) An ab plane filled by H-bonded chains 
running along [110]. (d) Next higher (or lower) ab plane filled 
by H-bonded chains running along [110]. In panels (b-d), all 
CI atoms, H atoms bonded to oxygens, and O atoms double- 
bonded to carbons have been omitted for clarity. 
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allel to the b axis because the two-fold screw symmetry 
cancels out the a and c components of the polarization 
of each chain. 

Interestingly, the proton is found to be almost midway 
between the nitrogen and oxygen. The observed N-H 
bond length of 1.41 A is much larger than the typical 
value of 1.03 A in proton-transfered ionic H-bondsJ^ It 
seems that the proton potential between oxygen and ni- 
trogen has only a single minimum, suggesting that Phz- 
H2ca might not be an order-disorder ferroelectric (FE), 
but rather a displacive-type one, i.e., where the paraelec- 
tric (PE) phase has a polar instability associated with a 
soft phonon mode. However, a dynamic proton fluctu- 
ation (i.e., a rapid back-and- forth motion of the proton 
between two minima located closer to N or O) has been 
suggested^ based on nuclear spin relaxation-time mea- 
surements using the 35 CI nuclear quadrupole resonance 
(NQR). The activation-energy barrier for proton transfer 
is estimated to be 68meV from the Arrhenius temper- 
ature dependence of the fluctuations £ Further support 
for this picture comes from a second ferroelectric phase 
(FE-II) that appears upon further cooling below 136 K 
(after passing through an incommensurate phase at 136- 
146 K). In FE-II, the proton is found to be completely 
transferred to the Phz nitrogen atomi The observed N- 
Hbond length of 1.12 A is now consistent with (although 
a bit longer than) a typical N-H + bond length (« 1.03 A) 
in other organic molecular salts 

In this work, we investigate the structures, energet- 
ics, spontaneous electric polarizations, lattice instabil- 
ities, and energy barriers for proton transfer by using 
first-principles calculations. There are two difficulties to 
a theoretical treatment of this important class of com- 
pounds. First, in order to predict stable crystal struc- 
tures and their properties, it is critical to include van 
der Waals interactions, which are important for inter- 
molecular interactions but are missing in conventional 
exchange-correlation (XC) functionals. Second, because 
of the small mass of the proton, the proton quantum fluc- 
tuations are large enough to significantly affect the rel- 
ative stabilities of different structures. To address these 
issues, we use the recently developed van der Waals den- 
sity functional (vdW-DF2) of Ref. [H and include the 
zero-point energy at the harmonic level. The results are 
compared with experiments as well as with calculations 
carried out using the semilocal functional of Perdew, 
Burke, and Ernzerhof (PBE)^ one of the most success- 
ful generalized-gradient functionals. We show that the 
proton transfer from H^ca to Phz lowers the energy by 
102 meV/unit-cell, and has an energy barrier of 8meV. 
No unstable phonon is associated with the proton trans- 
fer. 

The paper is organized as follows. In Sec. [II] we de- 
scribe the details of the computational methods used in 
the calculations. Then, in Sec. IIII| we present the re- 
sults of our calculations of lattice constants, energies, and 
phonons of polar and nonpolar structures, and of the en- 
ergy barrier for proton transfer. Finally, we summarize 



the work in Sec. [IVl 



II. METHODS 

We use the plane-wave pseudopotential method^ as 
implemented in the quantum-espresso package and 
Troullier-Martins norm-conserving pseudopotentialsJ^ 
We adopt an kinetic-energy cutoff of 80 Ry and a 2 x 6 x 1 
fc-point mesh for the Brillouin-zone sampling. All calcu- 
lations are done fully self-consistentl y 19 i 2Q using Soler's 
efficient algorithm 21 to treat the vdW-DF2 exchange- 
correlation energy functional. Atomic positions and lat- 
tice parameters are fully optimized until the residual 
forces and stresses are smaller than 7.7 me V/A and 0.5 
kbar, respectively. The Berry-phase technique^ 2 - is used 
to calculate the polarization of the structures. 

The transition path for the proton transfer and its 
energy profile are determined using the climbing-image 
nudged elastic band metho d 23 ! 24 with 7 images. All 
atoms are relaxed during the process, with the lattice 
parameters fixed to the experimental equilibrium values 
for the PE phase. 

The zero-point energy corrections are included at the 
harmonic level using the computed phonon frequen- 
cies. For the phonon calculation we use the experi- 
mental lattice parameters, the linear-response density- 
functional perturbation theor y 25 ! 26 for PBE, and the 
finite- difference method for vdW-DF2. The acoustic sum 
rule is imposed on the force constants. In order to make 
the phonon frequencies well converged up to 1 cm -1 , very 
tight convergence thresholds are used for the phonon 
calculation: tr2_ph and conv.thr are set to 10 -18 and 
10 -10 , respectively, two orders of magnitude smaller than 
typical values for inorganic solids. 

The FE and PE states have computed DFT energy 
gaps of 0.5 eV and 1.2 eV, respectively. (Since DFT 
tends to underestimate gaps, the true gaps are presum- 
ably larger.) All bands are fully occupied and there are 
no unpaired electrons. Thus, Phz-H2ca is clearly a band 
insulator, compatible with our choice of methods. 



III. RESULTS AND DISCUSSION 

A. Lattice constants 

First we validate our computational approach by cal- 
culating the lattice constants and comparing them with 
known experimental values as shown in Table HI The 
vdW-DF2 shows excellent agreement with experiments. 
For the PE phase, the deviations are 0.00 A (0%), 0.06 A 
(1%), and —0.14 A (—1%) along a, 6, and c, respectively. 
The relative deviations with respect to the experiments 
are given in the parenthesis. The deviations for the FE 
phase are similar: 0.05 A (0%), 0.09 A (2%), and -0.27 A 
(-2%). 
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FIG. 2. (Color online) Possible proton transfer (PT) processes in the unit cell, as indicated by large unfilled arrows, (a) 
Structure FE1, in which a H2ca donates a proton to a nitrogen in a neighboring Phz, leaving begind an oxygen lone pair on the 
H2ca and creating an electric dipole along the PT direction. The direction of polarization P is shown at the bottom left-hand 
corner as a thin arrow labeled as P; here the PT occurs along [110] so that P has both a and b components, (b) Structure 
FE2b, in which another PT occurs along [110] in the other H-bonded chain; P points along b by symmetry, (c) Structure FE2a, 
in which the second PT occurs along [110] instead; P points along a by symmetry, (d) Structure PE2, a doubly protonated 
paraelectric structure, (e)-(g) Enlarged view of the PT process. 



On the other hand, PBE, which is one of the most 
successful semi-local functionals, overestimates the PE 
lattice constants by -0.70 A (-6%), 1.08 A (28%) and 
—0.26 A (—2%) along a, 6, and c axis, respectively [for 
FE, -0.33A (-3%), 0.90A (24%) and -0.17 A (-1%)]. 
Except for the c lattice constant for the PE phase, all 
lattice constants are poorly reproduced. 

This comparison between vdW-DF2 predictions and 
experiments confirms that the vdW-DF2 functional is ca- 
pable of capturing all three important interactions (co- 
valent bonds, H-bonds, and van der Waals interactions) 
with good fidelity in this organic crystal. 



TABLE I. Comparison of experimental lattice constants with 
those calculated using PBE and vdW-DF2. 



Structure 



Axis 



Lattice constant (A) 







Expt. a 


PBE 


vdW-DF2 




a 


12.42 


11.73 


12.42 


Paraelectric 


b 


3.85 


4.93 


3.90 




c 


16.98 


16.72 


16.84 




a 


12.42 


12.10 


12.47 


Ferroelectric 


b 


3.79 


4.69 


3.88 




c 


16.91 


16.74 


16.64 



a Neutron diffraction experiments of Ref. measured at 300 K 
and 160 K for paraelectric and ferroelectric phase respectively. 



B. Proton-transferred structures 

In the PE phase, there is no proton transfer and all 
molecules are neutral. There are two H2ca molecules per 
unit cell, and each H^ca has two hydrogen bonds. The 
protons in those four hydrogen bonds can be transferred 
to neighboring Phz molecules. In this section, we con- 
sider all the possible proton transfer configurations con- 
sistent with the primitive-cell periodicity, starting from 
the simplest and working toward more complex ones. 

The simplest single-proton transfer is shown in 
Fig. Efa), where the H2ca at the center of the unit cell 
donates a proton to a neighboring Phz (to the left in 
this figure). The large arrow represents the direction of 
the proton transfer. We denote the resulting structure 
as FE1, indicating that it is ferroelectric and only one 
proton has transferred in the unit cell. 

This structural change breaks the centrosymmetry, 
making it ferroelectric, and lowers the energy by 137 meV 
per unit cell with respect to the paraelectric phase, which 
we denote henceforth as PE0. The relative energies of 
these and other structures (to be discussed shortly) are 
illustrated in Fig. [3l (All structures and energies in this 
section are calculated using vdW-DF2. A comparison 
with PBE will be given later.) The electric polarization 
P of the FE1 structure points approximately in the [110] 
direction as is shown by the thin arrow at the bottom 
left-hand corner of the figure. 

The [Hica] - molecule that already donated one pro- 
ton also has the possibility to donate a second one to the 
other neighboring Phz in the chain, as shown in Fig.[2jd). 
We denote this structure as PE2, where the '2' indicates 
that two protons are transferred in the unit cell. This 
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FIG. 3. (Color online) Energies of competitive proton- 
transferred structures presented relative to the FE2b ground- 
state energy, which is taken as zero. 



double protonation of a Phz restores the centrosymme- 
try, so that PE2 has no polarization. However this second 
proton transfer increases the energy by 309 meV, i.e., the 
energy of PE2 is 172 meV higher than that of PE. There- 
fore, in searching for the ground state, we do not consider 
other configurations based on double-protonation chains, 
i.e., those involving three (FE3) or four (FE4) transferred 
protons per cell. 

Two possibilities then remain, in which the other H^ca 
molecule at the corner site in the figure also donates a 
proton, as shown in Figs. [2jb-c). We found that the ex- 
perimental low-temperature ferroelectric structure that 
we already discussed in detail in the previous section, 
shown in Fig. Efb), has the lowest energy of all possi- 
ble configurations. We now denote it as FE26, indicating 
that two protons are transferred and the net polarization 
is parallel to the '6' axis. FE2b is more stable than PE0 
by 318 me V. This energy reduction is 44meV larger than 
twice the PEO-to-FEl reduction of 137 meV mentioned 
above. Thus, 44meV can be taken as an estimate of the 
interchain coupling strength along c. 

The polarization of FE2b is calculated to be 
4.5/iC/cm 2 , in a reasonable agreement with the mea- 
sured value of ~2 jiC / cm 2 . Also our calculated N-H bond 
length of 1.06 A is consistent with the experimental value 
of 1.12 A measured at low temperature in phase FE-II, 
and with a typical N-H + bond length of ~ 1.03 A in other 
organic molecular salts Thus, the structure of the fer- 
roelectric phase agrees also well with experiments. 

In the other case, in which the proton is transferred 
to the left Phz as in Fig. [2fc), only the a component of 
the polarization remains by the symmetry. Following the 
same notational scheme introduced above, we denote this 
as FE2a. Its energy is only 11 meV/unit-cell higher than 
that of the ground-state FE26 structure. 
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FIG. 4. (Color online) Energies relative to PE0 (left), as 
calculated with PBE and vdW-DF2, for the transition-state 
structure leading from PE0 to FE1 (middle) and for FE1 
(right). The lattice constants are fixed to experimental values 
for this comparison, and zero-point energies are included at 
the harmonic level. 



C. Energy barrier for proton transfer and lattice 
instability 

Here we investigate the proton transfer process from 
PE0 to FE1 in detail. The transition path for the proton 
transfer is calculated using the climbing-image nudged 
elastic band method j 23 i 24 A close-up view of the initial, 
transition, and final states of this proton transfer process 
are shown in Figs. [2fe)-(g). Along the transition path, 
the two molecules get closer to each other by up to 0.18 A. 
They then retreat again as the proton completes its trans- 
fer, but not completely; the final proton-transferred pair 
is closer by 0.05 A than the initial neutral one. 

It is well known that local or semi-local functionals 
such as PBE underestimate proton-transfer barriers j 27 i 28 
For example, for the intramolecular proton transfer in 
the malonaldehyde molecule, the energy barrier in PBE 
is only aquater of the known accurate value of 177 meV 
(Ref. l27|) calculated by the coupled-cluster method 
with single, double, and perturbative triple exciations 
[CCSD(T)]. On the other hand, the vdW-DF2 barrier of 
183 meV in that case is in excellent agreement with the 
accurate value, corroborating the validity of our method. 
For the Phz-E^ca crystal, we found a similar result; our 
computed PBE barrier of 44 meV from PE0 to FE1 is less 
than half of the vdW-DF2 value (105 meV). Furthermore, 
after including the zero-point corrections, the vdW-DF2 
barrier drops to 8 meV, and the transition becomes bar- 
rierless in PBE as is shown in Fig. |4j 

The energetics of proton transfer in vdW-DF2 is 
in good agreement with experimental observation of 
the thermally activated proton fluctuation in the high- 
temperature FE-I phase and the proton-transferred 
structure in the low-temperature FE-II phase. Also the 
apparent single- well proton potential in FE-I could be un- 
derstood by dynamic proton transfer in a small-barrier 
double-well potential. 
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Nevertheless, there could be other pathways not cap- 
tured by our initial NEB path, or an unstable phonon 
mode that could trigger the proton displacement toward 
the FE-I phase. However, our calculation shows that the 
zone-center phonons of the PEO structure are all stable. 
We do not find any unstable mode associated with the 
proton displacement toward the FE-I phase, the signa- 
ture of a displacive-type ferroelectric. 

Both the proton transfer energetics and the zone-center 
phonons are in good agreement with experimental obser- 
vations and support the picture of an order-disorder, as 
opposed to a displacive, FE transition. 

IV. CONCLUSION 

By using first-principles density-functional theory, we 
have studied the structure and energetics of a ferroelec- 
tric molecular crystal of phenazine and chloranilic acid, 
and have analyzed the energy barrier for proton transfer 
and the stability of lattice vibrational modes. We have 
shown that the inclusion of van der Waals interactions 
is crucial for a proper description of this molecular crys- 



tal, and that a recently developed vdW-DF2 functional 
reproduces the structures of the PE and FE phases in 
good agreement with experiment. We have found that 
the zone-center phonons of the PE state are all stable, 
and the proton transfer — the key structural change that 
leads to the ferroelectric structure — has an energy bar- 
rier of 8meV. The signature of a displacive-type ferro- 
electric, i.e., lattice instability in the PE phase, has not 
been found. Accordingly we propose that Phz-H 2 ca is 
an order-disorder FE. Our analysis of the stability of the 
lattice vibrational modes and the energy barrier for pro- 
ton transfer supports the possibility that the apparent 
single-well proton potential in the FE-I phase is an av- 
erage effect arising from dynamical proton transfer in an 
asymmetric double- well potential. 
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